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•  We  study  the  connection  between  fractional  calculus  and  subordination  processes. 

•  A  fractional  trajectory  is  equivalent  to  the  ensemble  average  trajectory  resulting  from  subordination. 

•  The  internal  friction  effects  associated  to  fractional  derivatives  can  alternatively  be  explained  by  decorrelation. 

•  We  use  subordination  theory  to  study  the  relaxation  toward  the  equilibrium  state  induced  by  the  fractional  derivatives. 

•  As  an  example  we  apply  the  theory  to  the  fractional  Lotka-Volterra  ecological  model. 


ARTICLE  INFO 


ABSTRACT 


Article  history: 

Received  17  April  2013 

Received  in  revised  form  21  June  2013 

Available  online  27  July  2013 


Keywords: 

Fractional  calculus 
Subordination 
Mittag-Leffler 
Nonlinear  operator 
Decorrelation 
Relaxation  to  equilibrium 


The  fundamental  connection  between  fractional  calculus  and  subordination  processes  is 
explored  and  affords  a  physical  interpretation  of  a  fractional  trajectory,  that  being  an 
average  over  an  ensemble  of  stochastic  trajectories.  Heretofore  what  has  been  interpreted 
as  intrinsic  friction,  a  form  of  non-Markovian  dissipation  that  automatically  arises  from 
adopting  the  fractional  calculus,  is  shown  to  be  a  manifestation  of  decorrelations  between 
trajectories.  We  apply  the  general  theory  developed  herein  to  the  Lotka-Volterra  ecological 
model,  providing  new  insight  into  the  final  equilibrium  state.  The  relaxation  time  to  achieve 
this  state  is  also  considered. 
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1.  Introduction 

In  the  last  few  years  there  has  been  an  increasing  interest  in  the  adoption  of  fractional  calculus  in  time  to  study  a  wide 
set  of  issues,  ranging  from  economics  [1-3]  to  ecology  [4,5]  to  allometry  [6]  to  bioengineering  [7[.  Fractional  calculus  more 
readily  caught  the  attention  of  engineers  since  it  allows  for  the  simulation  of  and  control  over  a  more  diverse  range  of 
complex  system  behaviors  [8],  The  primary  motivation  for  the  use  of  time  fractional  derivatives  has  been  that  they  generate 
solution  trajectories  with  memory,  and  therefore  give  a  more  realistic  representation  to  many  complex  processes  found  in 
nature.  This  memory  effect  has  produced  what  is  usually  referred  to  as  an  intrinsic  damping  in  the  system  [9-11],  pointed  out 
by  many  authors  to  be  equivalent  to  assuming  that  the  process  under  study  is  in  contact  with  a  dissipative  environment.  A 
fundamental  example  of  this  so-called  intrinsic  friction  effect  happens  with  the  fractional  simple  harmonic  oscillator,  where 
replacing  the  ordinary  time  derivatives  with  fractional  derivatives  leads  to  damped  oscillations,  even  without  considering  an 
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external  friction  force  [  10,11  ].  Perhaps  this  interpretation  of  the  fractional  harmonic  oscillator  stems  from  the  success  of  the 
fractional  calculus  description  of  viscoelastic  material  properties,  where  memory  and  internal  friction  are  plausible  [12-15], 

The  goal  of  the  present  article  is  to  suggest  the  importance  of  another  interpretation,  that  of  decorrelation.  We  show  that 
a  single  fractional  trajectory  can  be  represented  as  the  average  over  an  infinite  number  of  ordinary  differential  trajectories, 
each  trajectory  keeping  an  internal  time  using  stochastic  clocks  that  operate  independently  of  one  another  [16],  When  this  is 
the  case  the  dissipation  effects  noticed  by  the  majority  of  investigators  in  this  field  are  shown  herein  to  actually  be  a  form  of 
decorrelation.  To  support  our  theoretical  arguments  we  provide  an  ecological  example  in  an  application  of  the  theory  to  the 
Lotka-Volterra  model.  It  has  to  be  stressed,  however,  that  the  new  perspective  is  not  limited  to  ecology,  and  it  is  expected 
to  transform  the  interpretations  currently  adopted  in  many  other  fields  of  research,  as  discussed  in  the  concluding  section. 

Historically,  the  correct  physical  interpretation  of  fractional  calculus  has  been  elusive.  Perhaps  the  most  well  known 
example  of  a  dynamical  process  leading  to  a  fractional  calculus  description  lies  in  the  anomalous  diffusion  processes 
governed  by  the  continuous  time  random  walk  (CTRW)  theory  [17-19].  From  a  generalized  master  equation  approach, 
it  has  been  shown  that  an  ensemble  of  particles  undergoing  a  CTRW  in  an  external  field  can  be  represented  by  a  time- 
fractional  Fokker-Planck  equation  [20].  The  CTRW  is  a  specific  example  of  a  much  more  general  theory,  usually  referred 
to  as  subordination  [21-23]  because,  in  a  certain  sense,  subordination  is  a  diffusion  process  along  a  prespecified  path.  We 
emphasize  that  the  theory  presented  here  is  general  and  goes  beyond  diffusion  to  include  any  trajectory  resulting  from  the 
integration  of  fractional  differential  equations  in  time. 

In  Section  2  we  provide  a  general  demonstration  of  the  new  perspective  on  fractional  calculus  using  an  operator  pro¬ 
cedure  and  the  Mittag-Leffler  function.  Section  3  is  an  application  of  the  theory  to  the  Lotka-Volterra  system,  including  a 
comparison  between  the  fractional  Lotka-Volterra  trajectory  and  the  average  trajectory  resulting  from  the  subordination 
process.  In  addition  we  discuss  the  time  required  to  relax  to  the  equilibrium  state  induced  by  the  adoption  of  fractional  cal¬ 
culus.  Finally  Section  4  is  devoted  to  discussing  the  conceptual  consequences  that  the  new  fractional  derivative  perspective 
may  have  in  other  fields  of  investigation. 

2.  General  method 

In  this  section  we  demonstrate  the  equivalence  between  a  fractional  trajectory  that  is  the  solution  of  a  Caputo  fractional 
differential  equation,  and  the  ensemble  average  trajectory  that  results  from  a  subordination  process.  We  consider  only 
fractional  derivatives  of  the  Caputo  type,  noting  that  a  Riemann-Liouville  approach  would  be  equivalent  as  long  as  the 
initial  conditions  are  properly  specified. 


2.1.  Fractional  trajectory 


First  let  us  consider  the  fractional  differential  equation 

~j~V(t)  =  OV(t),  (1) 

where  0  <  a  <  1  and  0  is  an  operator,  either  linear  or  nonlinear,  acting  on  the  vector  V(t).  With  this  general  notation 
we  can  describe  any  dynamical  system  of  interest.  To  solve  Eq.  (1)  we  use  the  Laplace  transform  of  a  Caputo  fractional 
derivative  [24], 

=  s“V(s)  —  s“_1V(0),  (2) 

to  find  that 

V(S)  =  - ^V(O).  (3) 

s  -  Os1  “ 

This  is  the  Laplace  transform  of  a  Mittag-Leffler  function,  and  moving  back  to  the  time  domain  yields  the  fractional  trajectory 
V(t)  =  Ea(Ota)V(0).  (4) 

Note  that  the  Mittag-Leffler  function  is  defined  by 


df‘ 


:V(t);s 


r(ak  +  1)  ’ 

where  r  denotes  the  gamma  function  [19].  With  this  series  representation  it  can  be  seen  that  there  is  no  ambiguity  in  the 
action  of  the  operator  in  the  argument  of  the  Mittag-Leffler  function,  since  it  is  similar  to  the  exponential  solution  of  an 
ordinary  derivative  in  that  it  requires  the  calculation  of  a  series  of  operator  moments  of  the  form  0"V(0). 

Recall  the  main  idea  is  that  the  fractional  trajectory  described  by  Eq.  (4),  which  is  the  solution  to  the  fractional  differential 
equation  of  Eq.  (1),  is  equivalent  to  a  subordination  process.  To  establish  this  equivalence  we  will  now  move  to  the  process 
of  subordination. 


Ea(z)  =  J2 
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2.2.  Subordination  process 

Subordination  implies  the  existence  of  two  different  notions  of  time  [23].  One  is  the  operational  time  r,  which  is  the 
internal  time  of  a  single  unit,  with  a  unit  generating  the  ordinary  dynamics  of  a  non-fractional  system.  The  other  notion  is 
experimental  time  t;  the  time  as  measured  by  the  clock  of  an  external  observer.  Typically  in  the  operational  time  frame  the 
temporal  behavior  of  a  unit  is  regular  and  evolves  exactly  according  to  the  ticks  of  a  clock.  Therefore  it  is  assumed  that  the 
system’s  trajectory  in  operational  time  is  well  defined  and  given  by  V(t),  which  is  the  solution  to  the  ordinary  differential 
equation 


d 

— V  =  OV. 
dr 


(6) 


The  discrete  solution  to  Eq.  (6)  can  be  found  recursively  through 
V(r  +  At)  =  (1  +  OAr)V(r), 


(7) 


provided  that  the  integration  time  step  is  sufficiently  small  compared  to  the  eigenvalues  of  the  operator  0,  a  condition  which 
we  informally  state  as 


(8) 


OAx  «;  l. 


Moving  from  the  initial  condition  of  the  system  and  replacing  the  continuous  operational  time  r  with  its  discrete  counterpart 
n  we  can  write 


V(n)  =  (1  +  OAr)nV(0), 


(9) 


with  the  concise  notation 


V(n)  =  V(nA  r). 


(10) 


In  operational  time  the  units’  behavior  appears  ordinary,  but  to  an  experimenter  observing  the  units  their  temporal 
behavior  appears  erratic,  evolving  in  time  then  abruptly  freezing  in  different  states  for  extended  time  periods.  Because  of 
the  random  nature  of  the  experimental  time  evolution  of  the  units,  the  subordination  process  involves  an  ensemble  average 
over  many  units  each  evolving  according  to  its  own  internal  clock,  independent  of  one  another.  Making  an  ensemble  average 
over  a  large  number  of  units  results  in  a  smooth  average  trajectory,  which  is  equivalent  to  the  fractional  trajectory. 

To  find  the  average  behavior  we  can  move  from  the  operational  time  solution  to  the  experimental  time  solution  using 
the  probability  of  the  occurrence  of  an  event 


(11) 


n=0 


Actually  here  we  should  write  V(t)  as  (V(t)),  but  for  convenience  we  drop  the  brackets,  it  being  evident  that  the  trajectory 
resulting  from  the  subordination  process  inherently  involves  an  ensemble  average. 

To  interpret  the  physical  meaning  of  Eq.  ( 11 ),  consider  each  tick  of  the  internal  clock  n  of  a  unit  measured  in  experimental 
time  as  an  event.  Since  the  observation  is  made  in  experimental  time,  the  time  intervals  between  events  are  random.  We 
assume  that  the  waiting  times  between  consecutive  events  are  identically  distributed  independent  random  variables.  The 
integral  in  Eq.  (11)  is  then  built  up  according  to  renewal  theory  [25],  After  the  nth  event  the  unit  changes  from  state  V(n) 
to  V(n  +  1),  where  it  remains  until  the  action  of  the  next  event.  The  sum  over  n  takes  into  account  the  possibility  that  any 
number  of  events  could  have  occurred  prior  to  an  observation  at  experimental  time  t.  The  events  occur  probabilistically  with 
a  waiting-time  probability  distribution  function  (pdf)  i/r(t)  and  survival  probability  '/'(t).  The  waiting-time  pdf  is  related  to 
the  survival  probability  through 


(12) 


Taking  advantage  of  the  renewal  nature  of  the  events,  the  waiting-time  pdf  for  the  nth  event  in  a  sequence  is  connected  to 
the  previous  event  by 


(13) 


To  find  an  analytical  expression  for  the  behavior  in  experimental  time  it  is  convenient  to  study  the  Laplace  transform  of 


Eq.(ll), 


OO 


V(s)  =  $(s)£v(n)[£(s)]n, 


n=0 


(14) 
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where  the  relation  i/rn(s)  =  [0(s)]"  was  used  according  to  Eq.  (13).  With  the  discrete  time  solution  in  operational  time, 
Eq.  (9),  this  can  be  written  as 


V(s)  =  £(s)  ^(1  +  OAz)nV( 0)  [ir(s)f  . 


Performing  the  sum  and  noting  the  relationship  given  by  Eq.  (12)  we  find 
1  —  i/r(s)  1 


V(s)  = 

s  1  —  (1  +  OAzjifjs) 
This  can  be  expressed  in  the  typical  form 
1 


V(0). 


V(s)  = 


where 


4>(s) 


s  —  OAt<P(s) 

sty(s) 


V(0), 


1  -  \fr(s) 

is  the  Montroll-Weiss  memory  kernel  [17], 


(15) 


(16) 


(17) 


(18) 


2.3.  Establishing  an  equivalence  between  the  subordination  process  and  the  fractional  trajectory 


The  subordination  process  can  be  connected  to  the  fractional  trajectory  when  the  survival  probability  for  the  events  is 
given  by  the  Mittag-Leffler  function 


#(t)  =  Ea(-Ata),  (19) 

where  the  asymptotic  power-law  region  of  the  Mittag-Leffler  function  scales  with  the  same  index  a  as  the  fractional  deriva¬ 
tive  in  Eq.  (1).  In  order  for  the  Mittag-Leffler  function  to  be  compatible  with  a  survival  probability  the  index  must  be  in  the 
interval  0  <  a  <  1,  which  sets  the  interesting  condition  that  the  mean  waiting  time  between  consecutive  events  is  infinite 
since 

/»oo  /»oo 

(t)  =  I  r0(T)dr  =  /  £„(— /4r“)dr  (20) 

Jo  Jo 

diverges  due  to  the  fact  that  asymptotically  the  Mittag-Leffler  function  is  an  inverse  power-law.  This  inverse  power-law 
or  scale-free  condition  corresponds  to  nonstationary  processes  for  0  <  a  <  1  and  leads  to  a  form  of  ergodicity  break¬ 
down  [26],  Ergodicity  breakdown  has  interesting  statistical  consequences  at  the  level  of  single  trajectories  [27,28]  and  has 
been  observed  experimentally  in  anomalous  diffusion  processes  in  living  cells  [29,30], 

With  the  choice  of  survival  probability  given  by  the  Mittag-Leffler  function  the  Montroll-Weiss  memory  kernel  is  exactly 

0(s)=As'~°‘,  (21) 

and  Eq.  (17)  reduces  to 


V(s)  = 


1 

- ; — V(0). 

s  —  OArAs^~a 


Consequently,  under  the  condition 


AtA  =  1. 


(22) 

(23) 


Eq.  (22)  becomes  identical  to  the  Laplace  transform  of  the  fractional  derivative  and  the  solution  is  again  the  Mittag-Leffler 
function  given  in  Eq.  (4).  The  result  is  that  a  fractional  trajectory  can  be  represented  as  the  ensemble  average  trajectory  of 
a  subordination  process.  This  result  is  exact  when  the  survival  probability  for  events  is  a  Mittag-Leffler  function  and  the 
conditions  of  Eqs.  (8)  and  (23)  are  satisfied.  The  requirement  of  a  survival  probability  that  is  a  Mittag-Leffler  function  can 
be  relaxed  to  generic  survival  probabilities  that  share  the  asymptotical  inverse  power-law  behavior  of  the  Mittag-Leffler 
function,  by  way  of  a  stochastic  central  limit  theorem  [31], 


2.4.  Eigenfunction  expansion 

Let  us  examine  the  condition  defined  by  Eq.  (8)  more  carefully.  The  generic  system  described  by  Eq.  (6)  involves  a 
spectrum  of  eigenvalues  yj  of  the  operator  0.  The  eigenfunctions  associated  with  this  problem  are  given  by 

<Pj  (V)  =  (V,  Xj) 


(24) 
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where  Xj  are  the  eigenfunctions  of  the  adjoint  operator  Ch 

OfXj  =  -YjXj  (25) 

normalized  such  that  the  inner  product 

(vj,  Xk)  =  hk  (26) 

where  v,  is  an  eigenvector  of  0.  The  dynamics  of  the  eigenvectors  in  operational  time  are  given  by 

^  Xj)  =  (OV,  Xj)  =  (V,  OtXj.)  =  -rj  (V,  Xj)  =  -YA-  (27) 

which  has  the  exponential  solution 


4>i (V,  r)  =  e-hT0.(vo).  (28) 

Consequently,  the  state  vector  has  the  eigenfunction  decomposition,  as  long  as  the  operator  0  has  a  complete  spectrum, 

V(z)  =  J2  (V,  Xj)  v,  =  ^  Vj <t>j  (V,  r) .  (29) 

j  j 

We  can  select  a  discrete  set  of  points  along  the  continuous  curve  defined  by  the  eigenfunctions  separated  by  a  time 
distance  At  to  obtain 

0j(V,  nAr)  =  „),  (30) 

where  we  have  replaced  the  continuous  subordination  time  r  with  the  discrete  subordination  time  n.  The  subordination  of 
this  single  eigenfunction  following  the  previous  argument  results  in  the  exact  expression 


0,(V,s) 


1 

s  +  (1  —  e~YiAr)<t >(s) 


<Pj(V  o). 


Only  when  the  time  step  is  sufficiently  small  that  the  approximation 


(31) 


e~yJAT  sa  1  -  yjAt  (32) 

is  valid  over  the  entire  spectrum  of  eigenvalues  of  the  operator  0  can  we  write  Eq.  (22).  When  this  condition  is  satisfied  we 
recover  the  Laplace  transform  of  a  Mittag-Leffler  function  that  shares  the  properties  generated  by  the  fractional  trajectory 
in  Eq.  (3). 

This  digression  provides  the  interpretation  of  Eq.  (8).  If  the  discrete  time  step  is  properly  selected  then  the  fractional 
trajectory  will  be  exactly  equivalent  to  the  average  trajectory  resulting  from  the  subordination  process.  For  nonlinear 
systems  there  are  infinitely  many  eigenvalues  to  consider  when  choosing  the  time  step  At,  making  the  study  of  fractional 
trajectories  as  compared  to  subordination  a  non-trivial  problem  that  is  discussed  in  detail  in  the  following  section. 


3.  Ecological  application 

As  an  ecological  application  of  the  theory  presented  in  Section  2,  we  consider  the  Lotka-Volterra  equations  that  were 
originally  introduced  to  describe  the  population  dynamics  of  a  predator-prey  system.  The  fractional  system  we  propose 
strictly  serves  a  tutorial  purpose,  and  has  not  been  connected  to  any  real  ecological  observations  that  we  are  aware  of.  To 
construct  a  more  robust  fractional  system  with  a  stronger  connection  to  real  ecological  processes,  it  is  necessary  to  take  into 
account  a  proper  fractionalization  of  a  two-compartment  system,  see  for  example  [32,33]. 

This  section  is  divided  into  two  parts.  In  the  former  we  analyze  the  subordinated  Lotka-Volterra  system  and  compare 
it  to  the  Lotka-Volterra  fractional  trajectory.  In  the  latter  part  we  study  the  evolution  of  the  fractional  and  subordinated 
trajectories  toward  equilibrium. 


3.1.  A  comparison  between  the  fractional  trajectory  and  subordination 


The  fractional  Lotka-Volterra  system  is  described  by  [4,5] 
d“ 

— x  =  ax  —  bxy  (33) 

dt“ 

d“ 

— y  =  — cy  +  Sxy, 
d  t“ 


(34) 
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Fig.  1.  A  comparison  of  the  different  Lotka-Volterra  trajectories  in  the  predator-prey  phase  space.  The  ordinary  Lotka-Volterra  trajectory  in  operational 
time  (black)  is  a  cycle.  The  fractional  (blue)  and  average  subordinated  (green)  Lotka-Volterra  trajectories  spiral  in  to  the  fixed  point  (x  =  0.75,  y  =  0.5). 
An  ensemble  of  N  =  104  subordinated  trajectories  was  generated  to  find  the  smooth  average  trajectory  shown.  (For  interpretation  of  the  references  to 
color  in  this  figure  legend,  the  reader  is  referred  to  the  web  version  of  this  article.) 


where  x  represents  a  prey  population  and  y  the  predator  population.  According  to  the  general  arguments  of  the  previous 
section  we  can  define  the  operator  in  Eq.  (1)  as 

9  9 

0=  (ax-  bxy)  —  -  ( cy  -  5xy)  — ,  (35) 

3x  3y 

with  the  vector  V  having  the  components  x  and  y.  The  predictor-corrector  integration  method  is  adopted  to  solve  this 
system  of  fractional  differential  equations  to  generate  a  fractional  trajectory  when  the  fractional  index  is  a  =  0.9  [34].  For 
consistency  in  the  numerical  simulations  we  set  the  coefficients  to  the  values  a  =  0.001,  b  =  0.002,  c  =  0.003,  S  =  0.004, 
and  use  the  initial  conditions  (x0  =  0.4,  y0  =  0.6)  in  each  case.  This  fractional  trajectory  is  then  compared  to  the  ensemble 
average  trajectory  found  using  the  following  subordination  process. 

In  operational  time  the  ordinary  Lotka-Volterra  model  is  described  by  the  system  of  differential  equations 


— x  =  ax  —  bxy 
dr 

(36) 

d 

t -y  =  -cy  +  Sxy. 

dr 

(37) 

Using  the  predictor-corrector  method  with  a  =  1  we  numerically  integrate  the  system  of  differential  equations  to  find  the 
operational  time  trajectory  V(n)  =  (x(n),  y(n)),  where  we  have  set  the  operational  time  step  to  be  At  =  1.  This  ordinary 
Lotka-Volterra  trajectory  is  shown  in  Fig.  1.  We  then  go  through  the  subordination  process,  moving  from  operational 
time  to  experimental  time  by  generating  a  sequence  of  random  time  intervals  {Tj,  T2,  T3,  . . .},  from  a  waiting-time  pdf 
corresponding  to  a  Mittag-Leffler  function  survival  probability  [35],  The  parameters  of  the  Mittag-Leffler  function  were  set 
to  a  =  0.9  to  match  the  fractional  derivative,  and  A  =  1  to  fit  the  condition  given  by  Eq.  (23).  The  subordinated  trajectory 
begins  with  the  initial  condition  V(0)  =  (x0,  yo).  It  (the  trajectory)  stays  in  this  state  until  an  event  occurs,  at  experimental 
time  t  =  Ti  corresponding  to  the  first  random  waiting  time  drawn,  where  it  jumps  to  the  next  state  V(l).  Again  it  remains  in 
this  state  until  the  next  event  occurs  after  waiting  a  time  T2,  corresponding  to  an  overall  experimental  time  oft  =  T]  +  T2.  In 
this  way  a  single  subordinated  trajectory  is  constructed,  see  Fig.  2  for  example.  Applying  the  same  procedure,  an  ensemble 
of  subordinated  trajectories  is  created  and  from  it  the  average  over  subordinated  trajectories  is  found. 

A  comparison  between  the  ensemble  average  trajectory  of  the  subordination  process  and  the  fractional  trajectory  is 
shown  in  Fig.  1.  Both  the  fractional  trajectory  and  the  average  trajectory  resulting  from  the  subordination  process  spiral 
in  to  the  same  equilibrium  point,  however  the  two  trajectories  are  not  exactly  equivalent.  Here  we  are  forced  to  rest  on 
the  numerical  results  of  the  fractional  derivative  and  subordination  algorithms,  either  or  both  of  which  may  not  accurately 
represent  the  unknown  theoretical  trajectory.  If  we  make  the  unverified  assumption  that  the  fractional  derivative  algorithm 
produces  an  accurate  representation  of  the  fractional  Lotka-Volterra  system,  then  a  possible  source  of  inaccuracy  in  the 
subordination  algorithm  could  be  the  violation  of  the  condition  stated  in  Eq.  (8).  Due  to  the  nonlinear  nature  of  the 
Lotka-Volterra  system,  some  of  the  eigenvalues  of  the  operator  0  might  exceed  the  inverse  of  the  operational  time  step. 
The  large  eigenvalues  govern  the  short  time  behavior  of  the  fractional  and  subordinated  Lotka-Volterra  trajectories,  and  the 
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Fig.  2.  A  single  subordinated  Lotka-Volterra  prey  trajectory  observed  in  experimental  time  is  stochastic,  erratically  freezing  in  certain  population  levels 
for  times  that  can  be  on  the  same  order  as  the  overall  observation  time. 


difference  between  the  two  trajectories  may  be  due  to  an  inaccurate  representation  in  this  short  time  region.  In  fact,  both  the 
fractional  trajectory  and  the  average  trajectory  resulting  from  the  subordination  process  reach  the  same  equilibrium  point, 
an  indication  that  the  smallest  eigenvalues  of  the  operator  0,  determining  the  asymptotic  time  behavior,  fit  the  condition  of 
Eq.  (8).  For  nonlinear  systems  it  is  an  open  numerical  challenge  to  find  an  exact  equivalence  between  the  fractional  trajectory 
and  subordination  at  both  short  and  asymptotic  times,  whereas  for  linear  systems  an  exact  equivalence  between  the  two 
trajectories  can  be  found  for  all  times. 


3.2.  Decorrelation  and  equilibrium 


The  non-zero  fixed  point  for  the  ordinary  Lotka-Volterra  system  located  at  the  point 

(Xeq,yeq)  =  (^>  (38) 

can  be  found  by  setting  the  velocities  to  zero  in  Eqs.  (36)  and  (37).  This  fixed  point  remains  unchanged  when  moving  to  the 
fractional  Lotka-Volterra  system.  An  ordinary  Lotka-Volterra  trajectory  is  periodic,  making  one  cycle  after  a  period  T,  and 
the  fixed  point  can  alternatively  be  calculated  by  a  time  average  of  each  population  over  one  cycle  [36],  To  demonstrate  this 
calculation  for  xeq  we  go  back  to  Eq.  (37)  and  rewrite  it  as 

1  d 

y  =  -c  +  Sx  (39) 

y  dr 

which  has  the  solution 

In  =  —  ct  +  S  J  dtx(t).  (40) 

After  the  time  period  T  the  trajectory  returns  to  its  initial  condition,  y(T )  =  y0,  so  that  the  logarithm  on  the  left-hand  side 
of  Eq.  (40)  vanishes  identically  to  yield 


Xeq  — 


C 

S 


1 

T 


dtx(t). 


(41) 


These  calculations  can  be  repeated  to  findy,,,,  moving  from  Eq.  (36). 

This  time  average  can  be  transformed  into  an  ensemble  average  over  the  population  states.  To  make  this  transformation 
we  can  imagine  that  the  initial  condition  is  set  to  (xmjn,yeq),  where  xmin  is  the  minimum  value  that  the  x-population  can 
have  during  a  cycle.  The  trajectory  will  then  evolve  in  a  counter-clockwise  direction,  see  Fig.  1  for  example.  It  is  convenient 
to  divide  the  integral  in  Eq.  (41)  into  two  contributions 


f  \fo  dtX^t)  + jr  dtX^ 


(42) 
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Fig.  3.  The  probability  density  for  a  single  subordinated  Lotka-Volterra  trajectory  to  have  the  prey  level  x  at  experimental  times  q  =  103  (blue), 
t2  =  2  x  103  (red),  t3  =  3  x  103  (green),  and  t4  =  105  (black).  The  average  x-value  at  time  t4  is  xavg  =  0.7496,  indicating  that  p(x,  t4)  is  approximately  the 
equilibrium  distribution,  and  the  average  x-value  taken  from  the  equilibrium  distribution  is  equivalent  to  the  fixed  point  xeq  =  0.75.  N  =  106  units  were 
in  the  ensemble.  (For  interpretation  of  the  references  to  color  in  this  figure  legend,  the  reader  is  referred  to  the  web  version  of  this  article.) 


Here  V  is  the  time  corresponding  to  the  moment  when  the  lower  part  of  the  cycle,  where  y  <  yeq,  is  completed,  that 
is  x(T')  =  xmax,  the  maximum  value  that  the  x-population  can  have.  Now  we  make  a  change  of  variables  from  time  to 
population  space,  using  the  x-velocity,  vx  =  defined  in  Eq.  (36), 


HI 


X 


*(T') 

dx 

x(o>  vx(x,  y  ) 


+ 


rx(T) 

/  dx- 

Jx(T)  ! 


x(D  vx(x,y+) 


(43) 


The  x-velocity  is  positive  during  the  lower  part  of  the  cycle  where  y  <  yeq.  For  this  reason  we  label  the  x-velocity  in  this 
region  as  vx(x,  y~).  On  the  upper  part  of  the  cycle  where  y  >  yeq,  the  x-velocity  is  negative,  and  we  write  ux(x,y+)  to 
distinguish  this  case  from  the  lower  part  of  the  cycle.  Now  using  the  fact  that  x(0)  =  x(T)  we  can  combine  the  integrals  on 
the  right-hand  side  of  Eq.  (43)  to  get  the  result 


Xeq  — 


J  dxxpeq(x) 


where  the  equilibrium  pdf  depends  inversely  on  the  velocity  through 


(44) 


Peq  (x)  =\(  1  — - -H- — \  .  (45) 

T  \  vx(x,y  )  vx(x,  y+)  J 

Using  this  ensemble  averaging  approach,  the  equilibrium  point  reached  by  the  fractional  trajectory  can  be  understood 
intuitively  through  the  equivalent  subordination  process.  Imagine  an  ensemble  of  ordinary  Lotka-Volterra  systems  all 
identically  prepared  with  the  initial  condition  V(0)  =  (x0,  y0)  at  time  t  =  0.  Each  system  then  evolves  along  the  cycle,  but 
making  intermittent  steps  in  time  according  to  its  own  internal  clock.  The  probability  density  p(V,  f )  to  be  at  the  coordinate 
V  =  (x,  y)  at  time  t  spreads  from  the  initial  delta  function,  eventually  covering  all  possible  points  on  the  cycle,  as  some 
systems  step  ahead  while  others  get  stuck  for  extended  periods.  After  a  sufficiently  long  time  an  equilibrium  pdf  peq(V), 
that  is  proportional  to  the  inverse  of  the  x  and  y  velocities,  is  reached.  The  pdf  for  a  single  subordinated  trajectory  to  have  a 
prey  level  x  at  four  specific  experimental  times  is  shown  in  Fig.  3.  The  initial  delta  function  spreads  to  cover  the  entire  range 
of  possible  x-values,  and  tends  to  the  equilibrium  pdf  peq(x).  The  first  moment  of  the  prey  distribution  in  the  asymptotic 
time  limit  tends  toward  the  fixed  point  as  predicted  by  Eq.  (44).  An  equivalent  analysis  can  be  repeated  for  the  predator 
population. 

The  time  required  for  this  equilibrium  condition  to  be  reached  depends  on  the  index  a  of  the  fractional  derivative.  For 
a  =  1  there  is  no  decorrelation  between  the  identical  systems,  and  all  travel  together  along  the  cycle.  Therefore  we  can 
intuitively  expect  that  for  a  —>  1  the  time  to  reach  equilibrium  becomes  virtually  infinitely  extended  as  the  decorrelation 
among  the  systems  is  an  extremely  slow  process. 

Moving  from  a  discrete  operational  time  to  a  continuous  time  representation  [37],  Eq.  (11)  becomes 

POO 

V(t)=  /  dr  p(S,(t,  r)V(r) 

Jo 


(46) 
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where  p(S)(t,  r)  is  the  probability  that  the  operational  time  r  is  subordinated  to  the  experimental  time  t,  and  for  a  scaling 

pdf 


(47) 


The  time  derivative  of  the  trajectory  is  subsequently  given  by  the  approximate  expression 


(48) 


where  we  assume  V'(yt“)  can  be  replaced  by  V'(y)  fora  — »■  1. 

The  time  derivative  of  the  trajectory  tends  toward  zero  as  an  inverse  power-law,  but  for  a  close  to  one  it  may  take  a  very 
long  time  to  achieve  its  objective.  Theoretically,  a  small  region  in  phase  space  around  the  equilibrium  point  is  eventually 
reached  by  the  fractional  trajectory  for  any  a  <  1.  But  for  all  practical  purposes,  when  a  is  close  to  unity  the  decorrelation 
process  is  so  slow  that  the  equilibrium  condition  would  never  be  realized  in  numerical  calculations. 

4.  Concluding  remarks 

We  have  demonstrated  a  general  connection  between  fractional  calculus  and  the  process  of  subordination.  The  signifi¬ 
cance  of  this  connection  is  that  a  single  fractional  trajectory  obtained  as  the  solution  to  a  fractional  equation  of  motion  can 
be  interpreted  as  an  average  over  an  ensemble  of  stochastically  evolving  ordinary  trajectories.  This  interpretation  provides 
an  increased  physical  understanding  of  the  fractional  calculus.  Specifically,  by  adopting  a  decorrelation  perspective  rather 
than  the  notion  of  an  intrinsic  friction,  we  are  able  to  predict  the  time  it  takes  for  fractional  systems  to  reach  equilibrium 
using  subordination  theory. 

It  should  be  noted  that  accurate  numerical  integration  of  fractional  systems  is  an  area  of  active  research.  At  the  present 
time  there  is  no  conclusive  evidence  that  the  fractional  trajectory  resulting  from  the  predictor-corrector  integration 
accurately  represents  the  true  theoretical  fractional  trajectory.  We  leave  it  as  an  open  problem  to  determine  the  best 
numerical  approach  for  the  study  of  a  given  fractional  system,  which  would  involve  a  detailed  comparison  of  several  methods 
of  fractional  integration,  as  well  as  subordination.  Regardless  of  this  limitation,  the  work  presented  here  is  intended  to 
assist  in  improving  the  accuracy  of  the  current  numerical  methods.  For  example,  implementing  some  fractional  integration 
methods  may  result  in  a  fractional  trajectory  that  at  short  times  lies  outside  the  original  Lotka-Volterra  cycle.  This  effect 
appears  to  be  evidence  of  poor  accuracy.  It  cannot  be  realized  through  a  subordination  process,  where  decorrelation  always 
leads  to  a  trajectory  lying  completely  inside  the  original  cycle. 

The  new  insight  that  the  equivalent  subordination  process  brings  to  the  interpretation  of  fractional  trajectories  is  far- 
reaching  in  its  implications  for  complex  and  chaotic  systems.  It  is  well-known  that  memory  effects  can  suppress  chaos  [38], 
and  fractional  chaotic  systems  are  attracting  the  attention  of  many  researchers.  With  an  ensemble  average  perspective  the 
explanation  of  the  behavior  of  fractional  chaotic  systems  changes  dramatically.  Chaos  has  been  found  to  persist  in  fractional 
systems  when  a  remains  close  to  unity  [39].  On  the  basis  of  our  results  we  are  led  to  believe  this  may  be  a  consequence 
of  a  finite  observation  time.  If  the  observation  is  extended  infinitely,  we  expect  that  chaos  would  eventually  be  suppressed 
for  any  fractional  system  having  a  smaller  than  one.  Therefore,  adopting  the  new  perspective,  the  stability  analysis  of  these 
fractional  systems  can  be  significantly  improved  from  the  current  level  of  understanding.  This  is  moving  in  a  promising 
direction  as  the  connection  between  fractional  calculus  and  the  complex  processes  found  in  nature  continues  to  grow. 

Acknowledgments 

A.S.,  M.T.B.  and  P.G.  warmly  thank  Welch  and  ARO  for  their  support  through  Grant  Nos.  B-1577  and  W911NF-11-1-0478, 
respectively. 


References 


[1]  T.  Skovranek,  I.  Podlubny,  I.  Petras,  Modeling  of  the  national  economies  in  state-space:  a  fractional  calculus  approach,  Economic  Model.  29  (2012) 
1322-1327. 

[2]  F.  Mainardi,  M.  Raberto,  R.  Gorenflo,  E.  Scalas,  Fractional  calculus  and  continuous-time  finance  II:  the  waiting-time  distribution,  Physica  A  287  (2000) 
468-481. 

[3]  R.V.  Mendes,  A  fractional  calculus  interpretation  of  the  fractional  volatility  model,  Nonlinear  Dynam.  55  (2009)  395-399. 

[4]  S.  Das,  P.K.  Gupta,  A  mathematical  model  on  fractional  Lotka-Volterra  equations,  J.  Theoret.  Biol.  277  (2011)  1-6. 

[5]  E.  Ahmed,  A.M.A.  El-Sayed,  H.A.A.  El-Saka,  Equilibrium  points,  stability  and  numerical  solutions  of  fractional-order  predator-prey  and  rabies  models, 
J.  Math.  Anal.  Appl.  325  (2007)  542-553. 

[6]  D.  West,  B.J.  West,  On  allometry  relations,  Internat.  J.  Modern  Phys.  B  26  (2012)  1-57. 

[7]  R.L.  Magin,  Fractional  Calculus  in  Bioengineering,  Begell  House  Inc.,  Connecticut,  2006. 

[8]  D.  Cafagna,  Fractional  calculus:  a  mathematical  tool  from  the  past  for  present  engineers,  IEEE  Industrial  Electron.  Magazine  1  (2)  (2007)  35-40. 

[9]  V.E.  Tarasov,  Quantum  dissipation  from  power-law  memory,  Ann.  Physics  327  (2012)  1719-1729. 

[10]  A.  Tofighi,  The  intrinsic  damping  of  the  fractional  oscillator,  Physica  A  329  (2003)  29-34. 

[11]  B.N.  Narahari  Achar,  J.W.  Hanneken,  T.  Clarke,  Response  characteristics  of  a  fractional  oscillator,  Physica  A  309  (2002)  275-288. 


5672 


A  Svenkeson  et  al.  /  Physica  A  392  (2013)  5663-5672 


[12]  S.I.  Meshkov,  G.N.  Pachevskaya,  T.D.  Shermergor,  Internal  friction  described  with  the  aid  of  fractionally-exponential  kernels,  J.  Appl.  Mech.  Tech.  Phys. 
7(1966)63-65. 

[13]  F.  Mainardi,  Fractional  Calculus  and  Waves  in  Linear  Viscoelasticity,  Imperial  College  Press,  London,  2010. 

[14]  R.L.  Bagley,  P.J.  Torvik,  A  theoretical  basis  for  the  application  of  fractional  calculus  to  viscoelasticity,  J.  Rheol.  27  (1983)  201-210. 

[15]  Y.A.  Rossikhin,  M.V.  Shitikova,  Application  of  fractional  derivatives  to  the  analysis  of  damped  vibrations  of  viscoelastic  single  mass  systems,  Acta  Mech. 
120(1997) 109-125. 

[16]  A.A.  Stanislavsky,  Fractional  oscillator,  Phys.  Rev.  E  70  (2004)  051103. 

[17]  E.W.  Montroll,  G.H.  Weiss,  Random  walks  on  lattices.  II,  J.  Math.  Phys.  6  (1965)  167-181. 

[18]  R.  Metzler,  J.  Klafter,  The  random  walk’s  guide  to  anomalous  diffusion:  a  fractional  dynamics  approach,  Phys.  Rep.  339  (2000)  1-77. 

[  19]  B.J.  West,  M.  Bologna,  P.  Grigolini,  Physics  of  Fractal  Operators,  Springer,  New  York,  2003. 

[20]  R.  Metzler,  E.  Barkai,  J.  Klafter,  Deriving  fractional  Fokker-Planck  equations  from  a  generalised  master  equation,  Europhys.  Lett.  46  (1999)  431-436. 

[21]  H.C.  Fogedby,  Langevin  equations  for  continuous  time  Levy  flights,  Phys.  Rev.  E  50  (1994)  1657-1660. 

[22]  S.  Eule,  R.  Friedrich,  Subordinated  Langevin  equations  for  anomalous  diffusion  in  external  potentials— biasing  and  decoupled  external  forces,  EPL  86 
(2009)  30008. 

[23]  I.M.  Sokolov,  J.  Klafter,  From  diffusion  to  anomalous  diffusion:  a  century  after  Einstein’s  Brownian  motion,  Chaos  15  (2005)  026103. 

[24]  F.  Mainardi,  Y.  Luchko,  G.  Pagnini,  The  fundamental  solution  of  the  space-time  fractional  diffusion  equation,  Fract.  Calc.  Appl.  Anal.  4  (2001)  153-192. 

[25]  G.H.  Weiss,  Aspects  and  Applications  of  the  Random  Walk,  North-Holland,  Amsterdam,  1994. 

[26]  G.  Margolin,  E.  Barkai,  Nonergodicity  of  blinking  nanocrystals  and  other  Levy-walk  processes,  Phys.  Rev.  Lett.  94  (2005)  080601. 

[27]  S.  Burov,  R.  Metzler,  E.  Barkai,  Aging  and  nonergodicity  beyond  the  Khinchin  theorem,  PNAS  107  (2010)  13228-13233. 

[28]  Y.  He,  S.  Burov,  R.  Metzler,  E.  Barkai,  Random  time-scale  invariant  diffusion  and  transport  coefficients,  Phys.  Rev.  Lett.  101  (2008)  058101. 

[29]  J.-H.  Jeon,  V.  Tejedor,  S.  Burov,  E.  Barkai,  C.  Selhuber-Unkel,  K.  Berg-Sorensen,  L.  Oddershede,  R.  Metzler,  In  vivo  anomalous  diffusion  and  weak 
ergodicity  breaking  of  lipid  granules,  Phys.  Rev.  Lett.  106  (201 1)  048103. 

[30]  A.V.  Weigel,  B.  Simon,  M.M.  Tamkun,  D.  Krapf,  Ergodic  and  nonergodic  processes  coexist  in  the  plasma  membrane  as  observed  by  single-molecule 
tracking,  PNAS  108  (2011)  6438-6443. 

[31]  P.  Pramukkul,  A.  Svenkeson,  P.  Grigolini,  M.  Bologna,  B.J.  West,  Complexity  and  the  fractional  calculus,  Adv.  Math.  Phys.  2013  (2013)  498789. 

[32]  A.  Dokoumetzidis,  R.  Magin,  P.  Macheras,  Fractional  kinetics  in  multi-compartmental  systems,  J.  Pharmacokinetics  and  Pharmacodynamics  37  (2010) 
507-524. 

[33]  E.  Hanert,  Front  dynamics  in  a  two-species  competition  model  driven  by  Levy  flights,  J.  Theoret.  Biol.  300  (2012)  134-142. 

[34]  K.  Diethelm,  N.J.  Ford,  A.D.  Freed,  A  predictor-corrector  approach  for  the  numerical  solution  of  fractional  differential  equations,  Nonlinear  Dynam.  29 
(2002)  3-22. 

[35]  D.  Fulger,  E.  Scalas,  G.  Germano,  Monte  Carlo  simulation  of  uncoupled  continuous-time  random  walks  yielding  a  stochastic  solution  of  the  space-time 
fractional  diffusion  equation,  Phys.  Rev.  E  77  (2008)  021122. 

[36]  N.S.  Goel,  S.C.  Maitra,  E.W.  Montroll,  On  the  Volterra  and  other  nonlinear  models  of  interacting  populations,  Rev.  Modern  Phys.  43  (1971)  231-276. 

[37]  A.A.  Stanislavsky,  Probability  interpretation  of  the  integral  of  fractional  order,  Theoret.  and  Math.  Phys.  138  (2004)  418-431. 

[38]  A.A.  Stanislavsky,  Long-term  memory  contribution  as  applied  to  the  motion  of  discrete  dynamical  systems,  Chaos  16  (2006)  043105. 

[39]  I.  Grigorenko,  E.  Grigorenko,  Chaotic  dynamics  of  the  fractional  Lorenz  system,  Phys.  Rev.  Lett.  91  (2003)  034101. 


